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Abstract 

Non-equilibrium Molecular Dynamics (NEMD) calculations of the bulk viscosity of the triple 
point Lennard-Jones fluid are performed with the aim of investigating the origin of the observed 
disagreement between Green-Kubo estimates and previous NEMD data. We show that a careful 
application of the Doll's perturbation field, the dynamical NEMD method, the instantaneous form 
of the perturbation and the "subtraction technique" provides a NEMD estimate of the bulk viscosity 
at zero field in full agreement with the value obtained by the Green-Kubo formula. As previously 
reported for the shear viscosity, we find that the bulk viscosity exhibits a large linear regime with 
the field intensity which confirms the Lennard-Jones fluid as a genuine Newtonian fluid even at 
triple point. 
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I. INTRODUCTION 



The calculation of the hydrodynamics transport coefficients for model systems is a no- 
tl cea b le success of Molecule, D y ua mi cs (MD) j. The staudard m ethod to co m pu t e liuear 
transport coefficients by Molecular Dynamics simulations makes use of the Green-Kubo for- 
mulas . Based on the dissipation-fluctuation theorem, Green-Kubo formulas allow 
to compute linear transport coefficients from dynamical fluctuations of suitably defined mi- 
croscopic currents at equilibrium. The Green-Kubo methodology can be easily implemented 
in a simulation of the equilibrium state and all transport coefficients can be obtained in the 
same calculations. 

An alternative approach to the computation of transport coefficients is to mimic the 
appropriate non equilibrium state. This can be generally obtained by applying a suitable 
external force and measuring the response related to the corresponding transport coefficient. 
Non Equilibrium Molecular Dynamics (NEMD) has been developed along these lines already 
in the early eighties. It was soon realized that some paradigms of Equilibrium Molecular 
Dynamics (EMD) had to be relaxed to mimic non equilibrium processes. In particular the 
use of Periodic Boundary Conditions (PBC), a key ingredient of EMD to minimize finite size 
effects, is often incompatible with the non equilibrium state of interest. In many interesting 
cases the external field acts through the boundaries, for instance a thermal gradient or a 
velocity gradient, and the simulation of such system can require abandoning the use of PBC 
in favor of less convenient boundaries. This was for instance the case of a systems under 
the action of a thermal gradient a Poiseuille flow |7j. Non-periodic boundaries 

however require quite larger systems which had limited the early use of direct non equilibrium 
methods. To circumvent these limitations the so called "synthetic" NEMD algorithms have 
been developed and extensively used in the exploration of non equilibrium phenomena 8| . 
The general idea behind this class of algorithms is to replace the external force by an 
effective, PBC compatible, bulk field which, in the limit of vanishing intensity, excites the 
same response as the original external force. In this way the linear regime can, in principle, 
be explored without abandoning the use of PBC and therefore avoiding large finite size 
effects. In the case of fluid flows this technique requires the use of periodic but moving 
boundary conditions^, It should be noticed that, after restoring the use of PBC, the 
heat produced by the external bulk field must be removed by a "bulk" thermostatting 
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mechanism such as for instance a Nose- Hoover thermostat. We want to emphasize that the 
theoretical foundation of this class of algorithm is the Linear Response Theory and their use 
jeyond the linear regime is somewhat arbitrary. In this context the " subtraction technique" 
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, 112J is a very useful tool (at least for simple systems in which the response time is 
not longer than the typical Lyapunov time) to perform the vanishing perturbation limit. 
For almost all transport coefficients a good agreement between GK method and synthetic 



NEMD methods has been found 
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191 ] . The bulk viscosity makes a 



noticeable exception. In Fig. [T] we show the results of several computations of the bulk 
viscosity of a Lennard- Jones fluid close to the triple point. Most of these works adopted the 



Green-Kubo method 
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25]and found very similar results. Only two NEMD 



calculations have been performed so far 
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26j and they both provides values of the bulk 



viscosity 30%-50% higher than the Green-Kubo values. Note that finite size effects cannot 
explain the observed discrepancies. 

In the present work we reconsider the calculation of the bulk viscosity of the Lennard- 
Jones fluid close to the triple point and show that a careful application of the well know 
Doll's synthetic algorithm provides estimates of the bulk viscosity in full agreement with the 
Green-Kubo values. The paper is organized as follows. In section [Til we provide the necessary 
theoretical background by discussing both the Green-Kubo formula for the bulk viscosity 
coefficient and the dynamical NEMD approach we have adopted. Section II I II deals with 
details of the simulation such as the time dependence of the perturbation, the simulation 
box and the implementation of the dynamical approach to NEMD. In section HVl we collect 
our results and in section |V] some concluding remarks. In the Appendix we show that, 
in the linear regime, the synthetic perturbation for the viscosity is, as generally assumed, 
proportional to the velocity field produced. 



II. THEORETICAL BACKGROUND 



A. Hydrodynamics and microscopic identification of local fields 

The bulk viscosity is one of the transport coefficients introduced in hydrodynamics Q]. In 
this theory, the fluid is described by classical fields which, in the case of a simple neutral 
system, are the mass, momentum and energy density fields. Partial differential equations 
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FIG. 1: Review of previous results for the bulk viscosity of a simple Lennard- Jones fluid close to 
the triple point. The reported values are plotted as a function of the number of particles in the 
computations. 



derived from the continuity equation, supplemented by the so called "constitutive relations" 
and by the local equilibrium hypothesis, provide a closed theoretical framework for the evo- 
lution of these fields. The "constitutive relations" are linear relations between the external 
forces acting on the system and the excited flows, the coefficients being the transport coef- 
ficients specific for each material. The viscosity coefficients, namely the bulk viscosity, rj v , 
and the shear viscosity, 7] s , are defined by the Newton constitutive law 

£(f, t) = [p(f, t) - ( Vv - j^ s )V • v{r, t)]l - rj s [Vv(f, t) + (W(f, 0] (1) 

where P_ is the pressure tensor, v the velocity field and J the identity tensor, p represents the 
hydrostatic pressure, which, according to the local equilibrium hypothesis, can be expressed 
in terms of the mass density, m(f,t) = mn(f,t), and the energy density, e(r, t), by the 
equilibrium equation of state 

p{f, t) = P eq {mn(f, t), e(r, t)) (2) 

where n(f, t) is the local number density. If the velocity field reduces to the particular form 
v(f, t) = rf(t), with f(t) a yet unspecified function of time, the Newton law reduces to 

^Tr[E-pL] = -3 Vv f(t) (3) 

where the symbol Tr stands for the trace of the tensor. 
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According to Irving and Kirkwood [27(|, any macroscopic hydrodynamic field J(r,t) can 
be obtained from the statistical average, over the time dependent ensemble, p(T,t), of a 
microscopic observable J(T), where T = {ri,pi} (i = 1, N) is the phase-space point of the N 
particles system. In the specific case of the pressure tensor we havej^] 



N 



r — r,- 
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+ riFi 



(4) 



where is the internal force acting on particle i. The number and energy densities can also 
be expressed as 



N 



n(f, t) 



e(r, t) 



^Sif-Pi) p(T,t) 
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2$(f-fi) 
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2m 



p(r,t) 



(5) 
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where <p(r) represents the pair potential of our model system. The local hydrostatic pressure 
p(f, t) is the trace of the pressure tensor and, within the local equilibrium hypothesis, its 
fluctuations around the equilibrium value po{no, eo) can be expressed as 

dP f 



(e(f,t)-e ) + 
n on 



cq 



(n(r, t) - no) 



to 



(7) 



The bulk viscosity coefficient in terms of statistical averages is then obtained, for velocity 
fields of the form prescribed above, by replacing the fields in Eq.Q with the appropriate dy- 
namical ensemble averages. Furthermore, in order to ensure the validity of local equilibrium, 
it is required to take the "hydrodynamic limit" : 

Tr[P(k,uj)-p(k,uj)I] 



>1c 



— lim lim 



9/M 



(8) 



where J(k,uj) = h f f dr dtexp(ik ■ r) exp(iu}t)J(r,t) is the usual Fourier transform (in 
space and in time) of the field. 

B. Green-Kubo formula 



Equation (jSJ) expresses the bulk viscosity as an average on the non equilibrium distribution 
p(T,t). When the system is close to equilibrium (linear hydrodynamics) it is possible to 



rewrite it by the Green-Kubo formula 

Vv = PV 
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29|: 



dt < J(t)J(0) > 



(9) 



with 



J(t) = (V(t)- <V>)~y 



dP, 



eq 



de 



{H{t)- <H>) 



(10) 



where V, e are the volume and internal energy per unit volume of the system, respectively. 
V and ii. are the dynamical variables corresponding to the thermodynamic pressure and 
energy 
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V(t) = hm—Tr(E(k,t)) 
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J2 Tr 
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H(t) = hme(k,t) = J2 
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Eq. (TTU]) is the general expression of the current related to the bulk viscosity coefficients 
in the general case in which the energy fluctuates. If experiments are conducted in the 
microcanonic ensemble the energy fluctuations vanishes and the more familiar expression of 
the Green-Kubo formula is obtained 4|. 



C. "Doll's" perturbation and the 'dynamical approach' to Non Equilibrium 
Molecular Dynamics 

As described in the introduction, the alternative route to transport properties is to con- 
sider the system subjected to an external perturbation able to mimic the "thermodynamic 
force" which excites the appropriate nonequilibrium flux inside the system. In the present 
case the "thermodynamic force" is the macroscopic velocity gradient, Vv, while the corre- 
sponding flux is the deviation of the pressure from its local equilibrium value, |Tr[P] — p. 
The bulk external force to be used in such experiments is known as "Doll's" perturbation: 

N T 

H'{T,t) =^nPi ■■ (Vufat)) (13) 

i 

where V"u(r, t) is the required external field. This perturbation was proposed for the first 



time by Luttinger [30j and adopted in a Molecular Dynamics simulation by Hoover et al. [26, 
32J. In the Appendix we will show that, in the linear regime, the macroscopic velocity field 
v(f, t) induced by this perturbation coincides with imposed external constraint u(f, t) in the 
long wave-lenght limit 

~v(k = 0,u) = U(k = 0,u) (14) 



Once we are able to induce the required hydro-dynamic flux, we need a procedure to compute 
the average of the response on the non equilibrium ensemble. 



To this aim we can exploit the "Onsager-Kubo" relation[ll|, [12|, |31[. Calling S(t) the 
time evolution operator of the perturbed dynamics, the following relation holds for the non 
equilibrium average of the generic microscopic flux J 

<J>t = J J(T)p(T,t)dT = J J(T)S\t) Po (T)dT = 

= J S{t)J{T)p {T)dT = J J(t) Po (T)dT =< J(t) > (15) 

where S^(t) is the adjoint of S(t) and Po(T) is the ensemble distribution at the time t = 0. If 
the perturbation is switched on at time t = from an equilibrium state, po is the equilibrium 
distribution and the Onsager-Kubo relation (fl5|) allows us to compute the required average 
in a rigorous way Indeed, via standard equilibrium MD simulation, we can obtain a set of 
statistically independent configurations {1^} distributed according to po(T). Starting from 
those configurations we can follow the evolution of the system under the perturbed dynamics 
and obtain the required nonequilibrium average as the average of the evolved observable over 
the initial distribution according to the Onsager-Kubo relation. 

When a large perturbation is applied a thermostatting mechanism needs to be added 
to the equation of motion and the response can depend on it. Conversely for vanishingly 
small perturbations the standard form of linear response theory holds and the response 
depends only on the applied perturbation. In this limit however an additional numerical 
problem is encountered. The fluctuations of the microscopic variables are quite large and 
dominate the response in the limit of vanishing perturbations. In simple systems, where the 
response time is comparable to the Lyapunov time of the exponential divergence of nearby 
starting trajectories, the "subtraction technique" can be used to extract the signal out of 



the statistical noise 



11 



12J. If So(t) is the evolution operator representing the unperturbed 
dynamics with < S (t)J > = 0, the average values < S(t)J — S (t)J > and < S(t)J > 
are equal. On the other hand, their variances are very different: the thermal fluctuations 
of S(t)j and So(t)J are highly correlated and therefore largely cancel each others. This is 
true for times short enough. At large times the variance of the difference estimator becomes 
twice the variance of the simple estimator. 
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III. SIMULATION TECHNIQUE 



A. Impulsive external field 

As mentioned above, the current associated to the bulk viscosity to be used both in 
the Green-Kubo formula or in the NEMD experiments depends on the statistical ensemble 
chosen to conduct the experiment. In all previous equilibrium calculations, as well as in 
the present one, the microcanonical ensemble was chosen since the current reduces to the 
pressure tensor fluctuations without the need of evaluating the additional term related to 
the energy fluctuations. 

In NEMD, the Doll's perturbation for the bulk viscosity is a pure contraction or ex- 
pansion of the volume so that a constant perturbation in time, f(t) = e, will correspond 
to an exponential contraction/expansion of the volume (see next subsection), obviously an 
impractical way to extract the bulk viscosity coefficient. Alternatives forms of f(t) can be 
an oscillating function f(t) = esin(t) and an impulsive perturbation f(t) = eS(t), where 
e is the intensity of the field. The oscillating form, used in previous NEMD calculations, 
requires thermostatting the nonequilibrium trajectory in order to reach a steady state and 
the extended form of the current with the energy fluctuation term must be used. On the 
other hand, the impulsive perturbation acts on the system for an infinitesimal time, no ther- 
mostatting mechanism is necessary, and the dynamics after the impulse is the equilibrium 
dynamics for the isolated system. This fact greatly simplifies the NEMD experiment and 
the form of the flux, eq. ffTUl) . For each initial configuration Tj, the energy and the volume of 
the system change from their initial values 7io and Vq of the equilibrium system, to the time 
independent values 7i[ = fii{t > + ) and V = V(t > + ) (note that all replicas undergoes 
the same volume change). In order to obtain the appropriate flux in Eq. ffTUl) we have to 
calculate only the dynamical variable V(t) and its average asymptotic value: 

Poo = lim < V > t (16) 

where < ■ ■ ■ > t are averages over the non-equilibrium distribution at time t. Note that is 
not a properly defined pressure because the corresponding ensemble is not well defined since 
each member of the ensemble has a different energy. It is rather the asymptotic large time 
value of < V >t that needs to be subtracted to it in order to make the current integrable in 
time to provide the associated transport coefficient. 
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B. Periodic Boundary Conditions 



The explicit form of the Doll's perturbation field with a homogeneous velocity gradient 
Vit = f(t)I, is not compatible with the periodic boundary conditions of MD. However the 
periodicity of the system can be restored if we allow the box matrix H_ to evolve according 
to the external flow [33] 

K(t) = ■ M^) ( 17 ) 

Starting with a cubic cell of edge L a (0) = L Q ({a = x,y,z}), and applying a velocity field 
Vm = e8(t)l we get 

L a (t>0 + ) = L e*~L (l + e) (18) 
Similarly, the perturbation induces a discontinuity in the trajectory of the system 

fl = — + e8(t)fl 
m 

fi = Ft- eS(t)pi (19) 

which correspond to 

rl(0+) = r i (0-)e e ~r i (0-)(l + e) 

Pl(0 + ) = m-)e~ e ~ ft(0-)(l - e) (20) 

Therefore, the effect of the impulsive perturbation is to apply a homogeneous contrac- 
tion (expansion) of the position space and an homogenous expansion (contraction) in the 
momentum space of the system. Substituting eqs (120]) in the perturbed Hamiltonian 
7i(0 + ) = fio + ' P«(0 + ) an d using eq. (ECO), the variation of energy induced 

by the impulsive field is 

H(0 + ) - H(0-) = -V(0-)dV + C(e 2 ) 

where V is the volume of the system. Taking the ensemble average over the equilibrium 
distribution at t = 0~, we recover the first principle of thermodynamics 

AE = -pdV (21) 
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C. Bulk viscosity computation 



With the impulsive field of previous section, eq. (jHJ) for the bulk viscosity reduces to 



Vv 



lim lim 



lim 

e-+o 3e 



lTr[E~PL] 



Tr(Vv) 
<V>t ~Pc 



dt 



(22) 



Using the Onsager-Kubo relation to rewrite nonequilibrium ensemble averages in terms of 
equilibrium averages of evolved observables, and applying the subtraction technique de- 
scribed above, we obtain 



Vv 



lim 

e^o 3e 



< S(t)V - S (t)V > - [Poo - Po] dt 



(23) 



where S(t) and S (t) are the evolution operators for the perturbed and equilibrium dynamics 
respectively. Note that po = P eq (7io,no) is the pressure of the equilibrium microcanonical 
system, while as mentioned above, p^ is not a properly defined averaged pressure. 



D. Simulations scheme and numerical details 



As already mentioned, the operative procedure to compute eq. (1231 is to select a set of 
statistically uncorrelated equilibrium configurations of the system {Tj}, for each member of 
the set to generate both the equilibrium and the non equilibrium evolutions, and to compute 
the term [< S(t)V — So(t)V >o] up to a time t r , as the arithmetic mean over the set. The 
time t r is the typical time the individual non equilibrium system C; of energy 7ij(0 + ) takes to 
relax from the perturbed initial configuration to its equilibrium state. The additional offset 
term, involving p^, is the stationary value of that average beyond t r . In order to reduce 
the statistical noise on this term we have computed, for each individual nonequilibrium 
trajectory, the time average of the microscopic "pressure" between t r and 2 t r and we have 
estimated the offset as the arithmetic mean of those time averages over all nonequilibrium 
trajectories. 

The system considered in this work is a simple fluid of N = 864 particles interacting 
by the Lennard- Jones potential, in a thermodynamic state close to the triple point. In the 
following all quantities will be expressed in Lennard- Jones units: e = 1, a — 1, m = 1. The 
potential has been truncated at r cut = 2.5 and shifted to avoid discontinuity at r c . Moreover 
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FIG. 2: Simulations scheme: a long unperturbed trajectory of the system is used to sample initial 
conditions for segment of perturbed trajectories. Each integration of the perturbed equations 
of motion has been integrated for twice the estimated signal relaxation time. This is needed to 
evaluate the thermodynamic pressure of the segments. 

for r e [2.4 : 2.6] the potential has been replaced with a cubic polynomial in order to avoid 
discontinuities in the forces at the cutoff. The equations of motion have been integrated 
through a Velocity Verlet algorithm with an integration step h = 0.004436. 

The unperturbed trajectory (see figure [2]) has been integrated for 9.9 x 10 6 integration 
steps. No thermalizing device is added to the equilibrium dynamics so that a sampling of the 
microcanonical ensemble is obtained. This equilibrium trajectory was used to compute the 
bulk viscosity through the Green-Kubo formula. The time of saturation of the Green-Kubo 
integral, i.e. the decorrelation time of the pressure fluctuations at equilibrium, has been 
used as an estimate of the relaxation time t r . The set of initial equilibrium configurations 
are therefore selected as equilibrium configurations a time t r apart from each other along 
the equilibrium trajectory. In order to calculate the response of the system to the impulsive 
external field, S(t)V, the perturbed trajectories have been integrated for a time t r and 
extended to 2t r in order to evaluate the offset term (see figure [2]). 



IV. RESULTS 



The thermodynamic state of the equilibrium s yste m is re por ted in table [H This state is 



very close to that used in previous studies [lOj, |20j, |2l|, [22j, |23j, |24|, [25, 



26J. 
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TABLE I: Thermodynamic equilibrium state. The simulated system is in a state very close to that 
used in most of previous MD studies of the bulk viscosity. 

Density (r.u.) Energy/N (r.u.) Temp, (r.u.) Pressure (r.u.) 

0.8442 -4.112561(5) 0.72111(4) 0.8978(3) 



TABLE II: Fitting parameters for the tail of the Green-Kubo integral. 





fitting interval 


a(r.u.) 


6(r.u.) 


x 2 

ndf 


7] = at b 
rj = ae~ bt 


t e [0.22 : 2.0] 
t € [0.22 : 2.0] 


0.0160(2) 
0.1004(6) 


-1.09(1) 
1.83(1) 


le-05 
2e-06 



A. Green-Kubo results 

In Fig. [3] we report the Green-Kubo integrand R GK (t) = (3V < V(t)V(0) >o- After a 
rapid relaxation in about 0.22 time units, the curve exhibits a long time tail which vanishes 
only beyond t=2 (see the inset). Although the noise level on the time correlation function is 
quite small, the noise on its time integral, as obtained by a simple trapezoidal rule, results to 
be quite large because of the long tail. In order to get a smoothed signal we have attempted 
to replace the data beyond t=0.22 with two different analytic functions, fitted to the data in 
the time interval [0.22:2.0]. We have assumed a power law behavior R GK (t) ~ at~ b and an 
exponential behavior R GK (t) ~ a e~ bt (the power law behaviour is compatible with Hoover's 



hypothesis 261 ]: rj GK (u)) = a' + b'y/u, i.e. R GK (t) ~ t~?). Values of the fitting parameters 
are reported in table [Til while the data and the fitting functions are compared in figure HI 

From table [TTJ and Figure HJ we conclude that the exponential behaviour is a better 
representation of our data. Integration of the correlation function supplemented by our best 
exponential fit provides the behaviour in figure [5] and the following value for the viscosity 



Vv 



GK = 1.22 ± 0.03 



B. NEMD results 



In table [TTT1 we report the details and thermodynamic results from the performed nonequi- 
librium experiments. Note that positive e corresponds to nearly adiabatically expanded 
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FIG. 3: Green-Kubo integrand (see Eq.Q ). In the inset we show an enlargement of the tail. 
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FIG. 4: Tail of the Green-Kubo integrand R GK (t) and its fitting functions. Data are normalized by 
the initial value R GK (0). The fitting range is t £ [0.22 : 2]. In the left panel we show (in linear-log 
scale) the exponential fit, while in the right panel we show (in log-log scale) the power law fit. 
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FIG. 5: Green-Kubo integral: rf v (t) = Jq dsR GK (s). The horizontal plateau value is our Green- 
Kubo estimate of the bulk viscosity coefficient. 
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TABLE III: Thermodynamic properties of the system subjected to the external perturbation. As 
for the Poo, the values of the temperature T and of the energy variation AE are obtained as the 
arithmetic mean of the corresponding properties over all non equilibrium trajectories, see Section 
IIIIDI By comparing the data in the last two columns, we can verify the validity of the Eq. (l21l) , As 
explained in section IIIIBl this equation states the consistency of the '"Doll's" perturbation with 
the first principle of thermodynamics. 



e 


TV 


N/V 


T 


AE/N 


- P odV/N 


0.05 


9000 


0.7292 


0.5408(2) 


0.2191(1) 


0.2023(2) 


0.02 


9000 


0.7955 


0.6174(2) 


0.01049(6) 


0.02648(6) 


0.005 


13000 


0.8317 


0.6912(1) 


-0.01080(2) 


-0.02160(1) 


0.002 


16400 


0.8392 


0.70892(9) 


-0.005538(5) 


-0.006386(5) 


5xl0" 4 


16500 


0.8429 


0.71726(8) 


-0.001539(1) 


-0.001511(2) 


2xl0~ 4 


16500 


0.8437 


0.71909(9) 


-6.284(5) xl0~ 4 


-6.263(5) xl0~ 4 







0.8442 


0.72111(4) 






-2xl0~ 4 


12000 


0.8447 


0.7220(1) 


6.432(6) xlO" 4 


6.352(5) xl0~ 4 


-5xl0~ 4 


16500 


0.8455 


0.72335(9) 


0.001645(1) 


0.001660(1) 


-0.002 


16400 


0.8493 


0.7338(1) 


0.007222(5) 


0.006360(5) 


-0.005 


13000 


0.8570 


0.7533(1) 


0.02133(2) 


0.01585(1) 


-0.02 


9000 


0.8969 


0.86865(2) 


0.1603(1) 


0.06218(7) 


-0.05 


9000 


0.9846 


1.1894(5) 


0.9154(3) 


0.6539(5) 



systems and therefore to reduced temperatures (remember that the system remain isolated 
after the impulsive external perturbation at t=0). In the second last column we report the 
average variation of the total energy AE due to the impulsive perturbation. Comparing the 
data in the last two columns we can verify the validity of Eq. (|2"T|) in the limit of small e 
(|e| < 0.002). In order to analyze the response of the system, we need to separately discuss 
the various contribution to the integrand in eq. (1231) . For sake of clarity, let us define the 



14 



following quantities 



A(i) 

A°° 
R(t,e) 



S{t)V-S Q {t)T 



-3e 



Pc 



Po 



-3e 
A(t) - A°° 

dsR(s, e) 



(24) 

(25) 
(26) 
(27) 



The newtonian bulk viscosity is the zero field-infinite time limit of r] v (t, e). In table HVl and 
in figure El we report data for A(0) and A°° defined in Eqs. (1240 and ( 1251) . 

We note that the error on A°° grows when |e| decrease while the error on A(0) remains 
roughly constant. This is the effect of the subtraction technique that improves the signal to 
noise ratio for short t only while has no effect at large time where A°° has to be calculated. 
A less noisy estimate of A°° for |e| < 0.002 can be obtained by linear interpolation of the less 
noisy data at larger absolute values of the perturbation (see Fig. [6]). In Fig. [7] we show the 
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FIG. 6: Trend of the offset, A°°, defined in Eq. (l25p . as a function of the applied deformation e. 
For | e |< 0.002 the values are affected by a large error. An interpolation of the other data is used 
to reduce the noise in this region. 



values of the integrand in Eq.(J23J at t = + , i.e. R(0) = A(0) - A°°. At e = we display 
R GK (0). As predicted by Linear Response Theory, we observe that the NEMD response 
tends to the quadratic fluctuations of the pressure at equilibrium in the limit |e| — > 0. In 
Fig. [8] we show the estimates of the r] v (t,e) curves. We note that, as | e | decreases, the 
noise level at large time increases considerably. This signals again the limit of applicability 
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FIG. 7: Initial time response -R(0 + ) versus the perturbation. NEMD data tends to the correspond- 
ing Green Kubo value in the limit of small e. 



1.4 - 




^ 1.3 - 
% 1.2 - 




0.5 1 1.5 

time (r.u.) 



FIG. 8: Integral Tj v (t, e) of the response R(t, e), see Eqsj27l for some values of e. At large time and 
for small values of | e |, the subtraction technique is not able to reduce the noise in the signal. 
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of the subtraction technique. Similar to the Green-Kubo case of previous section, a less 
noisy estimator of the viscosity is obtained by integrating a response function in which the 
long time tail is replaced by an exponentially decaying behaviour fitted to the data at large 
time. Finally, in table HVl and in FigJHlwe report the values of r] v (e) = lim^oo rj v (t, e). As 
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FIG. 9: Values of t/„ calculated in the present work. In the inset we present a magnification of the 
small e region. A good agreement with the Green-Kubo result is obtained. 



expected the data in the small e region (|e| < 0.005) are in agreement with the Green-Kubo 
estimate of the viscosity. Although equilibrium and NEMD data are in agreement within 
error bars, rj v (e) data for |e| < 0.005 exhibit an unexpected small error bar and appear to 
be systematically below the Green Kubo prediction. This is probably a small bias of our 
extrapolation procedure. The large noise in the tail of the response function enforces us 
to perform the fit in a time interval considerably smaller than for larger perturbations (see 
table [TVT) . This might lead to underestimated errors (see FigJH]) and to an estimate of the 
asymptotic value slightly lower than the correct value. 

C. Non linear regime 

The results of previous section show the consistency between Green-Kubo and NEMD 
estimates of the bulk viscosity coefficient. However at variance with others coefficients, like 
for instance the shear viscosity, where a linear regime over several orders of magnitude of the 
intensity of the external perturbation is observed (up to roughly 0.05) |35J], in the present case 
the linear regime is apparently much reduced. This can be clearly seen in figure [7] where the 
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TABLE IV: In the second and in the third column we report the data for A(0) and A°°, respectively 
(see Eqs. (|24p and (|25p ) while, in the last column, the bulk viscosity for the corresponding intensity 
of the perturbation is shown. The time interval used for fitting the exponential function to the 
data is also reported in the second last column. 



f 


A(0) 


A°° 


range for the fit 




0.05 


25.132(5) 


13.757(6) 


0.44:2.66(500 pts) 


3.77(2) 


0.02 


32.884(5) 


21.382(7) 


0.44:2.66(500 pts) 


1.39(2) 


0.005 


37.778(4) 


25.61(2) 


0.44:2.66(500 pts) 


1.25(6) 


0.002 


38.863(4) 


26.53(4) 


0.44:2.22(400 pts) 


1.2(1) 


0.0005 


39.415(4) 


27.1(2) 


0.22:1.22(230 pts) 


1.15(6) 


0.0002 


39.527(4) 


26.6(4) 


0.22:1.22(230 pts) 


1.15(6) 


0(GK) 








1.22(3) 


-0.0002 


39.663(5) 


26.9(5) 


0.22:1.22(230 pts) 


1.13(6) 


-0.0005 


39.787(4) 


27.2(2) 


0.22:1.22(230 pts) 


1.15(6) 


-0.002 


40.351(4) 


27.80(4) 


0.44:2.22(400 pts) 


1.2(1) 


-0.005 


41.499(5) 


28.74(2) 


0.44:2.66(500 pts) 


1.18(5) 


-0.02 


47.888(6) 


34.116(7) 


0.44:2.66(500 pts) 


1.13(2) 


-0.05 


64.387(8) 


47.770(5) 


0.44:2.66(500 pts) 


1.15(2) 



NEMD results for R(0 + , e) matches the GK value with a finite slope suggesting that a linear 
expansion of the response function with the perturbation field is never justified. The same 
effect is seen for the viscosity in figure [9] (see the inset), although it is less pronounced and the 
much larger noise makes the observation less conclusive. In order to resolve this apparent 
paradox we must consider that the present perturbation, a contraction/expansion of the 
volume, excites the correct response and, at the same time, changes the thermodynamic 
state of the system. To correctly discuss the rheological non linear behavior of the fluid 
we should remove the pure thermodynamic contribution to the response. Let us consider 
the viscosity as a function of the perturbation, the internal energy and the volume of the 
system: r\ v = rj v (e,e,V). In the present case of impulsive perturbation we have ^(0 + ) = 
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TABLE V: Estimates of the thermodynamic derivatives in eqs. (|28p and (|29p at the considered 
state point as obtained by the central difference formula. 



dr/ v 

av ~ 


-0.0006(10) 


§| = -0.0024(8) 


8R(0) 
dV 


= -0.014(2) 


9 gP = 0.0021(9) 



H[l+3e+0(e 2 )], and£(0+) = E -p [V(0 + )-V }+O(e 2 ) = E -3ep V +O(e 2 ) where E ,V , 
and po represent the energy, volume and pressure of the equilibrium system respectively. The 
correct small e expansion for the viscosity is therefore 

e + 0(e 2 ) (28) 

A similar expansion holds for R(t), in particular for its value at t = + 

R(0 + ; e, e, V) = R(0 + ) e=0 + 

We have performed a series of EMD simulations at volumes and internal energies around 
the thermodynamic point studied and we have estimated the derivatives in equations (1281) . 
(|29|) by the central difference formula. In table |V] we report the estimated values of the 
derivatives. With those values the term in square brackets in eq. ( |29l) amounts to —50 ± 7 
to be compared with —59 ± 2, the value of the estimated slope of the response in NEMD 
data (see figure [TU]) . As for the viscosity itself, the value in the square brackets in eq (|2"8"|) is 
7±2 to be compared with 6±1 , the estimated slope of the viscosity in figure [9) When data 
for -R(0 + ) and r\ v are corrected by these thermodynamic terms a large linear regime appears 
as reported in figures [10] and [TlJ This behaviour confirms the simple LJ fluid as a true 
linear fluid in a large range of perturbations, the genuine rheological non linear behaviour 
appearing only beyond |e| ~ 0.02. 

V. CONCLUSIONS 

In the present paper we have reported Non-Equilibrium Molecular Dynamics calculations 
of the bulk viscosity of the Lennard- Jones fluid at triple point. Among the transport coeffi- 
cients of the simple fluid, the bulk viscosity was the only one for which NEMD results, from 



r] v (e,e,V) =r) v + 



P 



dE 



dR(0- 



v- 



d R(0 + ) \ 
~dE~) 



€ + 0{t 



(29) 



19 



15 
14.5 
14 
13.5 
c 13 
12.5 
12 



+ 
o 



-0.05 



12.8 
12.6 
12.4 
12.2 



corrected NEMD data 
NEMDdata 



-0.005-0.002 0.002 0.005 



-0.02 




£ 



0.02 



0.05 



corrected NEMD data 
NEMD data 



FIG. 10: Initial time response i2(0 + ).The NEMD data are corrected by the thermodynamic term 
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dr) v 



e=0 
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two independent previous studies, did not agree with the Green-Kubo estimates (see figure 
d]). Surprisingly, this unexpected failure of the Linear Response Theory remained unexplored 
for almost 25 years. In the present work we have resolved this apparent contradiction and 
found a full agreement between the NEMD and EMD estimates for the bulk viscosity. 

We have applied the Doll's perturbation field to excite the relevant flux in the system. 
At variance with the shear or elongational viscosity cases where the external perturbation 
is a superposition of a rotation and a deformation of the system at constant volume 3l|, |36) , 



the field needed to excite the flux related to the bulk viscosity coefficient is a compres- 
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sion/expansion of the volume at constant shape. Such a perturbation excites the desired 
flux but also changes the thermodynamic state of the system. As a consequence, the relevant 
flux depends on the conditions at which the non equilibrium experiment is conducted and, 



similarly, the Green-Kubo formula depends on the equilibrium ensemble used[29|. At equi- 
librium the microcanonical ensemble should be chosen to simplify the Green-Kubo analysis. 
As for the non equilibrium experiments, two different techniques can be applied. One can 
apply the stationary NEMD method in which the system is driver toward a stationary non 
equilibrium state by applying a periodic compression/expansion of the system. If the heat 
produced by the external field is removed by a thermostatting mechanism, the steady state 
can be maintained in time and the bulk viscosity can be estimated as the time average of the 
relevant flux divided by the perturbation strength. The other possible route is the dynam- 
ical NEMD in which a set of statistically uncorrelated replicas of the equilibrium system 
are subjected for t > to the perturbation field and the evolution of the ensemble can be 
followed in time under the perturbed dynamics. The dynamical method is superior to the 
stationary method because not only steady state informations can be obtained but also the 
transient behavior can be fully characterized. Moreover within the dynamical NEMD, the 
subtraction technique can be used to perform the zero field strength limit and to extract 
the value of the linear transport coefficient. Another advantage of the dynamical method is 
that one can apply an impulsive perturbation rather than a periodic one. Since the system 
is perturbed for a very short period of time (one step of our discrete dynamics) we do not 
need to introduce a thermostatting mechanism to perform a meaningful experiment. Using 
dynamical NEMD with the subtraction technique and an impulsive form of the perturbation, 
we were able to explore a large range of perturbation strength and carefully study the small 
field regime. We have found a perfect agreement between the NEMD and the Green-Kubo 
estimates of the bulk viscosity. These estimates are also in agreement with previous values 
obtained by the Green-Kubo formula for various system sizes, while they do not agree with 
previous NEMD studies conducted by the stationary NEMD method. Finally, by remov- 
ing the thermodynamic contribution to the viscosity, we show that the Lennard- Jones fluid 
at triple point exhibits a large linear regime, in agreement with the results for the shear 
viscosity |34J]. 
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APPENDIX: PERTURBATION ASSOCIATED TO VISCOUS FLUX: "DOLL'S" 



In this section we will compute the response in the velocity field to the "Doll's" pertur- 
bation. By a straightforward application of the linear response theory we will find the result 
in Eq. <HM . It proves that the "Doll's" perturbation is the perturbation producing every 
kind of viscous flux. 

Let us consider a system of iV particles of mass m with Hamiltonian of such a system 
will be: 

IV ^ N 

ff <0 ' = £!; + 5>(l ? «l) ( A -i) 

1=1 j^i 

where =f\ — fj and <p(r) is the pair potential. We add a perturbation term of the form: 

H {I) (T,t)= I g(f) : g(r,t)dr (A.2) 
Jv 

where V is the volume of the system, a(r) is a local dynamical variable and $(r, t) is the 
external field, dependent on time t and space r. We assume that $ is proportional to a 
small parameter e defining the magnitude of the perturbation. The linear response theory 
states that the effect of on a given observable O of the system is: 

poo r 

0(f, t) =< 6> -(3 dr ds < 6(0, 0)g(s, r)p >: g(f - s,t-r) + 0(e 2 ) (A.3) 
Jo Jv 

In the present case is given by EqJT^I therefore we set: a = Yli ?iV%&{fi — t) and 

We have to calculate the velocity field v for the system subjected to the perturbation. 
By following Irving-Kirkwood theory, we can identify this field by exploiting the equation: 

J m (r,t) = m(r,t)v(r,t) =< g(r) > t 

where we have introduced the flux of momentum J m and the corresponding dynamical 
variable g: 

N 

g(f) = mfiS(fi — r) (A. 4) 



Expanding the last expression in series of the small parameter e, we find: 

J m (r,t) M\f,t)+0(e 2 ) 



v(f, t) 



m(r, t) m(°) (f, t) + mW (f, t) + C(e 2 ) 
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mN 

where we have assumed that in the limit e — > 0: J m (f,t) — > and m(r,t) — > >n y-. Now, 
by applying the linear response theory Eq. flA.3j) . we find the following equation for the a 
component of the field v: 

v r°° r - 

v a (r, t) = J dr J ds < g a (0, 0)d>(s, r) > V ^u v {r + s,t - t) + 0(e 2 ) 

Finally, by performing an integration by parts we obtain: 

y poo p 

v a (r, t) = ^jjP J dr J ds < g a (0, Q)V^a^ u (s, r) > u v {r- s,t-r) + 0(e 2 ) 
which is a convolution in space and time. Its Fourier transform is 

v a (k, to) = cr a)3 (k, u)up(k, lo) (A. 5) 

where 

v r +oc r 

cr a p(k,uj) = —P J dr J dse li ~ k - g+ ^ < g a (0, 0)V M «^(s, r) > = 

= J dr < g a (0,0)(V^ lf3 )(k,r) > 

where we have performed the integral in ds in the second last equality. 

In order to evaluate a a /3(k,u), we have to calculate the Fourier transform (V^d^). In 
the following equation, we report the result of this calculation that will be demonstrated at 
the end of this paragraph. In the case with k << a, where "a" is the mean free path, we 
will find: 

N 

(V M d M/3 )(A : ,r) = -g p {k,t) - ik^ ■ • h){r$i) & e~ ihfi (A.6) 

i 

With the aid of the last equation, we can determine a a/ 3, defined in Eq. IA.5[ in the limit 
k -> 0: 

V f + °° ■ - 

\im a a/3 (k,uj) = —/3 / dr e tuJT < g a (0,0)g p (k = 0,r) > = 

V f +OD 

= -^jjP J dr iuoe WT < g a {0, 0)g^{k = 0, r) > (A.7) 
The ensemble average < g a (0, 0)g/3(k = 0, r) > can be calculated as follows: 
< 9a{r,0)gp{k^ 0,r) > =< g a (r,0) J dsgp{s,r) > = 
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N N 

=< ^2mr ia d(fi - r) ^ mr i/? (r) > 
* j 

where we have used the definition of g given by Eq JA.4[ Now, we note that the quantity 
Ylf m ^j/3( T ) corresponds to the total momentum of the system, therefore it is independent 
of the time. On the other hand, the overall average on the equilibrium ensemble have to be 
r independent as well. Therefore the following relation holds: 



1 f N N 

< g a (f, 0)gp(k = 0, r) > = — / df < ^ mr ia 5{fi - r) } j mr jf3 > 



N N N 

y<2_^ mr ^ mr iv >0= y < mr * a >0 al3 = m pv p 



J 

In the last equality we have also applied the equipartition theorem < mrf a > = K b T = i. 
By replacing the last result in Eq. (I A. 71) we finally find: 

/•+oo /"+oo 

Iimg(fc, u) = - dr iue iujT l = -iu / rfr e iuJT 9(r)l = 

JO J-oo 



where stands for the step function. By comparing this relation to Eq. ( 1A.5I) we obtain the 



final result v(0, uS) = u(0, lu) that proves the correspondence between the field u involved in 
the "Doll's" perturbation and the macroscopic velocity field v induced in the system by the 
perturbation. 

In order to complete the paragraph we have to demonstrate equation Eq. (1A.6I) . First of 
all, we calculate the derivative a as follows: 



^ • _ ( ■ d \ 

4(r, t) = ^(riPi + riPi)5(ri - r) + ( • t^t^ - *0 J 



t V j^i / 

where we have applied the relation -^rS(fi — r) = — V<5(r*j — r). 
Then, we calculate the Fourier transform on the space variable r: 



Q 



*) = J2 ~ fi E D + ^(f - ■ k) ) (A.8) 
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By focusing on the former equation term that involves the derivative of the potential, we 
find: 



N N 

EE 



d<t>(\ rg \) ( 
Of 



1 N N 

^EE 



i j^i 

. N N 



df 



+ To 



iEE 

1 N N 

^EE 



g0G gj I) . 

A g^(| gj I) 



e ik ' ?i ( fe 1 ^ - f 



N N 



EE 



r y + f,(ik ■ n, + 0(fc • f,; 2 ))] ®i 



By means of this result and under the hypothesis of small wave vector limit ka << 1, Eq. 
( 1A.8|) can be written as follows: 



v 



N 



a 



m\ ra I) 
df 



N 



+ ^{ik ■ fiWSi) e ik ' f 



(A.9) 



Eq. (1A.9P implies that the Fourier transform of the gradient of a fulfills the following 
equation: 

Jk-fj^ d0(| fjj 

13 df 



(V • a}(k,r) 



*• E 



mr,r,- e ! — 



iV 



ik -^{ik ■ fi){fiPi) e ik '^ 



This result correspond to the required Eq. (IA.6I) providing that the quantity in square 
brackets is equal to —g(k,t). This can easily verified as follows: 

N N 



d 



g(k,t) = — mfje lk ' ri = ymfi{ik ■ f) + mfj e 



8=1 



ik-fi 



1=1 
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AT / AT 

^ I mfiik-fi) - ^ 
i=l \ j^i 

. AT AT 



<9n 
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= ik ■ {mrif^j 



e ik-fi _ 1 ^2 ^ lJ { e ik -^ — e ik - ? 
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1=1 jyi 



)- 



E ifc • (m^/i) e** - * - - E E 



30(1 5j l) c< fc.j 



i=l i=l j'^i 

In the limit ka << 1 we finally get 

N N N 

i=l jV« 
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